Data Analysis for pathogenicity across variation in microbial virulence and host resistance

Load all required libraries

If one or more packages is not already installed for you, you can install it using install.packages("PackageName")

library(readr)
library(dplyr)
library(stringr)
library(survival)
library(ggplot2)
library(survminer)
library(wesanderson)
library(ggpubr)
library(coxme)
library(lme4)
library(corrplot)
library(RColorBrewer)
library(ggrepel)
library(reshape2)
library(ggcorrplot)
library(nlme)

Set your working directory and import dataset

Set your working directory to the folder where your dataset is present

The dataset can be downloaded from Main Data and Bacterial Load Data

setwd("~/Documents/GitHub/Pathogen-Potential")

Data <- read_csv("PathogenPotential_AG_04062023.csv")

LoadData <- read_csv("PathogenPotential_BacterialLoad_AG.csv")

Data processing and transformation to survival format

Entering Dpt host genotype information

Data <- Data %>%
  mutate(Genotype = recode(Genotype,
                           "1G" = "DptNull",
                           "1N" = "DptR",
                           "8N" = "DptS"))

Converting survival data to 0/1 format

#Create an empty dataframe Data_1 to store the survival information

Data_1 <- data.frame(Vial_ID = character(),
                     Genotype = character(),
                     Treatment = character(),
                     OD_group = character(),
                     Date = character(),
                     Day_Dead = numeric(),
                     censor = numeric())


#Convert data to 0/1 format

for(i in 1:length(Data$Vial_ID)){
  temp=data.frame(Data[i,])
  
  for(j in seq(10,22,2)){
    k = temp[1,j-2]-temp[1,j]
    if(k>0){
      for(kk in 1:k){
        Data_1[length(Data_1$Vial_ID)+1,]=c(temp[1,1],
                                            as.character(temp[1,2]),
                                            as.character(temp[1,3]),
                                            temp[1,5],
                                            as.character(temp[1,6]),
                                            j/2 - 4,
                                            1)
      }    
    } 
  }
  
  for(j in 8){
    k = temp[1,j]
    if(k>0){
      for(kk in 1:k){
        Data_1[length(Data_1$Vial_ID)+1,]=c(temp[1,1],
                                            as.character(temp[1,2]),
                                            as.character(temp[1,3]),
                                            temp[1,5],
                                            as.character(temp[1,6]),
                                            j/2 - 4,
                                            0)
      }    
    } 
  }
  
  for(j in 22){
    k = temp[1,j]
    if(k>0){
      for(kk in 1:k){
        Data_1[length(Data_1$Vial_ID)+1,]=c(temp[1,1],
                                            as.character(temp[1,2]),
                                            as.character(temp[1,3]),
                                            temp[1,5],
                                            as.character(temp[1,6]),
                                            j/2 - 4,
                                            0)
      }
    }
  }
}


#Converting Data_1 columns from character to numeric/factors

Data_1$Day_Dead = as.numeric(Data_1$Day_Dead)
Data_1$censor = as.numeric(Data_1$censor)

Data_1$Genotype <- as.factor(Data_1$Genotype)
Data_1$Treatment <- as.factor(Data_1$Treatment)
Data_1$Date <- as.factor(Data_1$Date)
Data_1$OD_group <- as.factor(Data_1$OD_group)

#Data_1 is the dataframe containing only survival data 

Calculating fraction symptomatic

Dead individuals on any day are assumed to be symptomatic

Data$Fs_Day1 <- (Data$Day0_Alive - Data$Day1_Climbing) / Data$Day0_Alive
Data$Fs_Day2 <- (Data$Day0_Alive - Data$Day2_Climbing) / Data$Day0_Alive
Data$Fs_Day3 <- (Data$Day0_Alive - Data$Day3_Climbing) / Data$Day0_Alive
Data$Fs_Day4 <- (Data$Day0_Alive - Data$Day4_Climbing) / Data$Day0_Alive
Data$Fs_Day5 <- (Data$Day0_Alive - Data$Day5_Climbing) / Data$Day0_Alive
Data$Fs_Day6 <- (Data$Day0_Alive - Data$Day6_Climbing) / Data$Day0_Alive
Data$Fs_Day7 <- (Data$Day0_Alive - Data$Day7_Climbing) / Data$Day0_Alive


#Creating a dataframe Fs for fraction symptomatic by transforming data
Fs <- melt(Data[,c(1:3,5,6,24:30)], id=c("Vial_ID",
                                         "Genotype",
                                         "Treatment",
                                         "OD_group",
                                         "Date"))
colnames(Fs)[7] <- "Fs"


Fs[c('F', 'Day')] <- str_split_fixed(Fs$variable, '_Day', 2)
Fs <- Fs[,c(1:5,7,9)]

Calculating Pathogen Potential

PP = (Fs/I) * (10^M)

Data$PP_Day1 <- (Data$Fs_Day1/exp(Data$OD_group)) * 
  (10 ^ ((Data$Day0_Alive - Data$Day1_Alive)/Data$Day0_Alive))

Data$PP_Day2 <- (Data$Fs_Day2/exp(Data$OD_group)) * 
  (10 ^ ((Data$Day0_Alive - Data$Day2_Alive)/Data$Day0_Alive))

Data$PP_Day3 <- (Data$Fs_Day3/exp(Data$OD_group)) * 
  (10 ^ ((Data$Day0_Alive - Data$Day3_Alive)/Data$Day0_Alive))

Data$PP_Day4 <- (Data$Fs_Day4/exp(Data$OD_group)) * 
  (10 ^ ((Data$Day0_Alive - Data$Day4_Alive)/Data$Day0_Alive))

Data$PP_Day5 <- (Data$Fs_Day5/exp(Data$OD_group)) * 
  (10 ^ ((Data$Day0_Alive - Data$Day5_Alive)/Data$Day0_Alive))

Data$PP_Day6 <- (Data$Fs_Day6/exp(Data$OD_group)) * 
  (10 ^ ((Data$Day0_Alive - Data$Day6_Alive)/Data$Day0_Alive))

Data$PP_Day7 <- (Data$Fs_Day7/exp(Data$OD_group)) * 
  (10 ^ ((Data$Day0_Alive - Data$Day7_Alive)/Data$Day0_Alive))

#Creating a dataframe PP for pathogen potential by transforming data
PP <- melt(Data[,c(1:3,5,6,31:37)], id=c("Vial_ID",
                                         "Genotype",
                                         "Treatment",
                                         "OD_group",
                                         "Date"))
colnames(PP)[7] <- "PP"

PP[c('p', 'Day')] <- str_split_fixed(PP$variable, '_Day', 2)

PP <- PP[,c(1:5,7,9)]

Calculating proportion survival

#Proportion survival on Day X
Data$Prop_Day0_Alive <- Data$Day0_Alive/Data$Day0_Alive
Data$Prop_Day1_Alive <- Data$Day1_Alive/Data$Day0_Alive
Data$Prop_Day2_Alive <- Data$Day2_Alive/Data$Day0_Alive
Data$Prop_Day3_Alive <- Data$Day3_Alive/Data$Day0_Alive
Data$Prop_Day4_Alive <- Data$Day4_Alive/Data$Day0_Alive
Data$Prop_Day5_Alive <- Data$Day5_Alive/Data$Day0_Alive
Data$Prop_Day6_Alive <- Data$Day6_Alive/Data$Day0_Alive
Data$Prop_Day7_Alive <- Data$Day7_Alive/Data$Day0_Alive

Calculating cumulative risk scores https://stat.ethz.ch/R-manual/R-devel/library/survival/html/predict.coxph.html

Data_Risk <- subset(Data_1, !Day_Dead==0)
cox_model <- coxph(Surv(Day_Dead, censor) ~ as.factor(Vial_ID) +
                     Genotype + 
                     Treatment + 
                     OD_group, 
                   data = Data_Risk)

# Compute the risk scores
risk_scores <- predict(cox_model, type = 'risk')

#Add it to a dataframe
Data_Risk$Risk_Score <- risk_scores

#group by vial id and day to assign cumulative scores

Data_1_2 <- Data_Risk %>%
  group_by(Vial_ID,Genotype,Treatment,OD_group) %>%
  summarise(Vial_ID=Vial_ID,
            Genotype=Genotype,
            Treatment=Treatment,
            OD_group=OD_group,
            Risk_Score=sum(Risk_Score)) %>%
  distinct()

#Add risk scores to the main dataframe
Data <- merge(Data, Data_1_2)

Calculating PP_T & Fs_T & Fs/T

T=median survival day

#PP_T Fs_T Fs/T

#Estimating mean survival day
Surv <- melt(Data[,c(1:4,6,39:45)], id=c("Vial_ID",
                                         "Genotype",
                                         "Treatment",
                                         "OD_group",
                                         "Date"))
colnames(Surv)[7] <- "Prop_alive"

Surv[c('p', 'Day')] <- str_split_fixed(Surv$variable, '_Alive', 2)
Surv[c('p', 'Day')] <- str_split_fixed(Surv$p, 'Prop_Day', 2)

Surv <- Surv[,c(1:5,7,9)]

Surv$Day_T <- 7

for(i in 1:length(Surv$Vial_ID)){
  temp=data.frame(Surv[i,])
  
  for(j in 6) {
    k=temp[1,j]
    
    if(k<=0.5) {
      Surv$Day_T[i] = Surv$Day[i]
    }
  }
}


Surv_1 <- Surv %>%
  group_by(Vial_ID,Genotype,Treatment,OD_group) %>%
  summarise(Vial_ID=Vial_ID,
            Genotype=Genotype,
            Treatment=Treatment,
            OD_group=OD_group,
            Day_T=min(Day_T)) %>%
  distinct()

Data <- merge(Data, Surv_1)

Data$Fs_T <- NA
Data$PP_T <- NA

#Estimating pathogen potential and fraction symptomatic on mean survival day T

for(i in 1:length(Data$Vial_ID)) {
  temp=data.frame(Data[i,])
  
  j=temp[1,47]
  
  if(j==1) {
    Data$Fs_T[i] = Data$Fs_Day1[i]
    Data$PP_T[i] = Data$PP_Day1[i]
  }
  
  if(j==2) {
    Data$Fs_T[i] = Data$Fs_Day2[i]
    Data$PP_T[i] = Data$PP_Day2[i]
  }
  
  if(j==3) {
    Data$Fs_T[i] = Data$Fs_Day3[i]
    Data$PP_T[i] = Data$PP_Day3[i]
  }
  
  if(j==4) {
    Data$Fs_T[i] = Data$Fs_Day4[i]
    Data$PP_T[i] = Data$PP_Day4[i]
  }
  
  if(j==5) {
    Data$Fs_T[i] = Data$Fs_Day5[i]
    Data$PP_T[i] = Data$PP_Day5[i]
  }
  
  if(j==6) {
    Data$Fs_T[i] = Data$Fs_Day6[i]
    Data$PP_T[i] = Data$PP_Day6[i]
  }
  
  if(j==7) {
    Data$Fs_T[i] = Data$Fs_Day7[i]
    Data$PP_T[i] = Data$PP_Day7[i]
  }
  
}

Data$`Fs_T/T` <- Data$Fs_T/as.numeric(Data$Day_T)

Main Figures

Figure 1B

Fig1b <- ggplot(Data, aes(x = as.factor(OD_group),
                 y = Prop_Day4_Alive,
                 pch = as.factor(OD_group),
                 col = Genotype,
                #fill = as.numeric(OD_group)
                )) +
  stat_summary(geom = "errorbar", 
               position = position_dodge(0.7), 
               width = 0.3,
               #aes(alpha = as.numeric(OD_group))
               ) +
  stat_summary(geom = "point", 
               position = position_dodge(0.7), 
               size = 6,
               #aes(alpha = as.numeric(OD_group))
               ) +
  facet_grid(~Treatment, scales = "free_x") +
  coord_cartesian(ylim=c(0,1)) +
  scale_x_discrete(
    labels = function(x) {
      if ("LB" %in% unique(Data$Treatment)) {
        # For LB, show only '0'
        ifelse(x == "0", "0", x)
      } else {
        # For other treatments, show the specific values
        ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
      }
    }
  ) +
  scale_shape_manual(values = c(4,1,2,0,5,7,8)) + 
  theme_classic() +
  labs(x = "OD",
       y = "Proportion alive on day 4",
       pch="OD") +
  scale_color_manual(values = c("red","blue","darkgreen")) +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 30),
        legend.title = element_text(size = 30),
        axis.title = element_text(size = 30),
        strip.background = element_blank(),
        strip.text = element_text(size = 30, face = "italic"),
        axis.text = element_text(size = 30))

Figure 1B - Interval plot of the proportion of flies alive on day 4 post infection for various host genotypes (homozygous serine dpt (wildtype, dptS69), homozygous arginine dpt (dptS69R), and dpt null (Δdpt)), and variable pathogen strains (P. burhodograneria Strain D, P. rettgeri, P. sneebia, P. alcalifaciens) with different loads of infective inoculum (OD).

Figure 2

Data_1$Treatment <- factor(Data_1$Treatment, levels = c("LB",
                                                        "P. alcalifaciens",
                                                        "P. burhodogranareia strD",
                                                        "P. rettgeri",
                                                        "P. sneebia"))

fit <- survfit(Surv(Day_Dead, censor) ~ OD_group + Treatment + Genotype,
               data = Data_1)

# Figure 2A
Survival <- ggplot(data = surv_summary(fit), 
                   aes(x = time, 
                       y = surv, 
                       group = interaction(Genotype,
                                           OD_group))) +
  facet_grid(~Treatment, scales = "free_x") +
  coord_cartesian(ylim=c(0,1)) +
  scale_x_discrete(
    labels = function(x) {
      if ("LB" %in% unique(Data$Treatment)) {
        # For LB, show only '0'
        ifelse(x == "0", "0", x)
      } else {
        # For other treatments, show the specific values
        ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
      }
    }
  ) +
  theme_classic() +
  labs(x = "Days post infection",
       y = "Survival") +
  geom_step(aes(color = Genotype),
            size = 1.5,
            alpha=0.4) +
  geom_point(aes(color = Genotype, 
                 shape = OD_group), 
             size = 4) +
  scale_shape_manual(values = c(4,1,2,0,5,7,8)) +
  coord_cartesian(xlim = c(0, 7)) +
  scale_color_manual(values = c("red",
                                "blue",
                                "darkgreen")) +
  labs(pch = "OD") +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 30),
        legend.title = element_text(size = 30),
        axis.title = element_text(size = 30),
        strip.background = element_blank(),
        strip.text = element_text(size = 30, face = "italic"),
        axis.text = element_text(size = 30))

Fs$Treatment <- factor(Fs$Treatment, levels = c("LB",
                                                "P. alcalifaciens",
                                                "P. burhodogranareia strD",
                                                "P. rettgeri",
                                                "P. sneebia"))
# Figure 2B
FracSick <- ggplot(subset(Fs,
                         Day==4), 
                aes(x = as.factor(OD_group), 
                    y = Fs, 
                    pch = as.factor(OD_group),
                    col = Genotype)) +
  stat_summary(geom = "errorbar", 
               position = position_dodge(0.7), 
               width = 0.3,
               #aes(alpha = as.numeric(OD_group))
               ) +
  stat_summary(geom = "point", 
               position = position_dodge(0.7), 
               size = 6,
               #aes(alpha = as.numeric(OD_group))
               ) +
  facet_grid(~Treatment, scales = "free_x") +
  coord_cartesian(ylim=c(0,1)) +
  scale_x_discrete(
    labels = function(x) {
      if ("LB" %in% unique(Data$Treatment)) {
        # For LB, show only '0'
        ifelse(x == "0", "0", x)
      } else {
        # For other treatments, show the specific values
        ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
      }
    }
  ) +
  scale_shape_manual(values = c(4,1,2,0,5,7,8)) + 
  theme_classic() +
  labs(x = "OD",
       y = "Fs on day 4",
       pch="OD") +
  scale_color_manual(values = c("red",
                                "blue",
                                "darkgreen")) +
  theme(legend.position = "top",
        legend.text = element_text(size = 30),
        legend.title = element_text(size = 30),
        axis.title = element_text(size = 30),
        strip.background = element_blank(),
        strip.text = element_text(size = 30, face = "italic"),
        axis.text = element_text(size = 30)) +
  guides(color=FALSE)

Data$Treatment <- factor(Data$Treatment, levels = c("LB",
                                                    "P. alcalifaciens",
                                                    "P. burhodogranareia strD",
                                                    "P. rettgeri",
                                                    "P. sneebia"))
# Figure 2C
Risk <- ggplot(Data, 
                aes(x = as.factor(OD_group), 
                    y = Risk_Score,
                    pch = as.factor(OD_group),
                    col = Genotype)) +
  stat_summary(geom = "errorbar", 
               position = position_dodge(0.7), 
               width = 0.3,
               #aes(alpha = as.numeric(OD_group))
               ) +
  stat_summary(geom = "point", 
               position = position_dodge(0.7), 
               size = 6,
               #aes(alpha = as.numeric(OD_group))
               ) +
  facet_grid(~Treatment, scales = "free_x") +
  scale_x_discrete(
    labels = function(x) {
      if ("LB" %in% unique(Data$Treatment)) {
        # For LB, show only '0'
        ifelse(x == "0", "0", x)
      } else {
        # For other treatments, show the specific values
        ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
      }
    }
  ) +
  scale_shape_manual(values = c(4,1,2,0,5,7,8)) + 
  theme_classic() +
  labs(x = "OD",
       y = "Fs on day 4",
       pch="OD") +
  scale_color_manual(values = c("red",
                                "blue",
                                "darkgreen")) +
  labs(x = "OD",
       y = "Risk Score",
       pch="OD") +
  theme(legend.position = "top",
        legend.text = element_text(size = 30),
        legend.title = element_text(size = 30),
        axis.title = element_text(size = 30),
        strip.background = element_blank(),
        strip.text = element_text(size = 30, face = "italic"),
        axis.text = element_text(size = 30)) +
  guides(color=FALSE)


Figure2 <- ggarrange(Survival, FracSick, Risk, 
                  nrow = 3, ncol = 1, 
                  common.legend = TRUE, 
                  labels = "AUTO",
                  font.label = list(size = 40))

Figure 2 - Host-based measures of pathogenicity:

• A. Survival curves,

• B. Fraction symptomatic on Day4, and

• C. Risk scores for infected flies of various host genotypes (homozygous serine dpt (wildtype, dptS69), homozygous arginine dpt (dptS69R), and dpt null (Δdpt)), and variable pathogen strains (P. burhodograneria Strain D, P. rettgeri, P. sneebia, P. alcalifaciens) with different loads of infective inoculum (OD).

Figure 3

# Figure 3A
PP_D4 <- ggplot(Data, 
                aes(x = as.factor(OD_group), 
                    y = PP_Day4,
                    pch = as.factor(OD_group),
                    col = Genotype)) +
  stat_summary(geom = "errorbar", 
               position = position_dodge(0.7), 
               width = 0.3,
               #aes(alpha = as.numeric(OD_group))
               ) +
  stat_summary(geom = "point", 
               position = position_dodge(0.7), 
               size = 6,
               #aes(alpha = as.numeric(OD_group))
               ) +
  facet_grid(~Treatment, scales = "free_x") +
  scale_x_discrete(
    labels = function(x) {
      if ("LB" %in% unique(Data$Treatment)) {
        # For LB, show only '0'
        ifelse(x == "0", "0", x)
      } else {
        # For other treatments, show the specific values
        ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
      }
    }
  ) +
  scale_shape_manual(values = c(4,1,2,0,5,7,8)) + 
  theme_classic() +
  labs(x = "OD",
       y = "Fs on day 4",
       pch="OD") +
  scale_color_manual(values = c("red",
                                "blue",
                                "darkgreen")) +
  labs(x = "OD",
       y = "PP on day 4",
       pch="OD") +
  theme(legend.position = "top",
        legend.text = element_text(size = 30),
        legend.title = element_text(size = 30),
        axis.title = element_text(size = 30),
        strip.background = element_blank(),
        strip.text = element_text(size = 30, face = "italic"),
        axis.text = element_text(size = 30))

# Figure 3B
PP_T <- ggplot(Data, 
                aes(x = as.factor(OD_group), 
                    y = PP_T,
                    pch = as.factor(OD_group),
                    col = Genotype)) +
  stat_summary(geom = "errorbar", 
               position = position_dodge(0.7), 
               width = 0.3,
               #aes(alpha = as.numeric(OD_group))
               ) +
  stat_summary(geom = "point", 
               position = position_dodge(0.7), 
               size = 6,
               #aes(alpha = as.numeric(OD_group))
               ) +
  facet_grid(~Treatment, scales = "free_x") +
  scale_x_discrete(
    labels = function(x) {
      if ("LB" %in% unique(Data$Treatment)) {
        # For LB, show only '0'
        ifelse(x == "0", "0", x)
      } else {
        # For other treatments, show the specific values
        ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
      }
    }
  ) +
  scale_shape_manual(values = c(4,1,2,0,5,7,8)) + 
  theme_classic() +
  labs(x = "OD",
       y = "Fs on day 4",
       pch="OD") +
  scale_color_manual(values = c("red",
                                "blue",
                                "darkgreen")) +
  labs(x = "OD",
       y = "PP_T",
       pch="OD") +
  theme(legend.position = "top",
        legend.text = element_text(size = 30),
        legend.title = element_text(size = 30),
        axis.title = element_text(size = 30),
        strip.background = element_blank(),
        strip.text = element_text(size = 30, face = "italic"),
        axis.text = element_text(size = 30))

# Figure 3C
Fs_T <- ggplot(Data, 
                aes(x = as.factor(OD_group), 
                    y = Fs_T,
                    pch = as.factor(OD_group),
                    col = Genotype)) +
  stat_summary(geom = "errorbar", 
               position = position_dodge(0.7), 
               width = 0.3,
               #aes(alpha = as.numeric(OD_group))
               ) +
  stat_summary(geom = "point", 
               position = position_dodge(0.7), 
               size = 6,
               #aes(alpha = as.numeric(OD_group))
               ) +
  facet_grid(~Treatment, scales = "free_x") +
  coord_cartesian(ylim=c(0,1)) +
  scale_x_discrete(
    labels = function(x) {
      if ("LB" %in% unique(Data$Treatment)) {
        # For LB, show only '0'
        ifelse(x == "0", "0", x)
      } else {
        # For other treatments, show the specific values
        ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
      }
    }
  ) +
  scale_shape_manual(values = c(4,1,2,0,5,7,8)) + 
  theme_classic() +
  labs(x = "OD",
       y = "Fs on day 4",
       pch="OD") +
  scale_color_manual(values = c("red",
                                "blue",
                                "darkgreen")) +
  labs(x = "OD",
       y = "Fs_T",
       pch="OD") +
  theme(legend.position = "top",
        legend.text = element_text(size = 30),
        legend.title = element_text(size = 30),
        axis.title = element_text(size = 30),
        strip.background = element_blank(),
        strip.text = element_text(size = 30, face = "italic"),
        axis.text = element_text(size = 30))

Figure3 <- ggarrange(PP_D4, PP_T, Fs_T, 
                  nrow = 3, ncol = 1, 
                  common.legend = TRUE, 
                  labels = "AUTO",
                  font.label = list(size = 40))

Figure 3

• A. Pathogen potential on Day 4 post-infection,
• B. Pathogen potential on median survival day (T), and
• C. Fraction symptomatic on median survival day (T) for infected flies of various host genotypes (homozygous serine dpt (wildtype, dptS69), homozygous arginine dpt (dptS69R), and dpt null (Δdpt)), and variable pathogen strains (P. burhodograneria Strain D, P. rettgeri, P. sneebia, P. alcalifaciens) with different loads of infective inoculum (OD).

Figure 4

Figure 4A

# Variables of interest
variables_of_interest <- c("Risk_Score", 
                           "Prop_Day4_Alive", 
                           "Fs_Day4", 
                           "PP_Day4", 
                           "Fs_T", 
                           "PP_T")

# Create an empty dataframe to store p-values matrices
Data_stats_1 <- data.frame(
  variable = character(),
  OD_group = character(),
  Genotype = character(),
  Treatment = character(),
  `OD_group*Genotype` = character(),
  `OD_group*Treatment` = character(),
  `Genotype*Treatment` = character(),
  `OD_group*Genotype*Treatment` = character(),
  stringsAsFactors = FALSE
)
TableS3 <- list()

# Loop through each variable of interest
for (variable in variables_of_interest) {
  # Fit an ANOVA model
  model <- aov(as.formula(paste(variable, "~ OD_group * Genotype * Treatment + Error(Vial_ID+Date)")), data = Data)
  
  # Extract ANOVA table
  anova_table <- as.data.frame(summary(model$Within)[[1]])
  
  TableS3[[variable]] <- anova_table
  
  # Create a new row to append
  new_row <- data.frame(
    variable = variable,
    OD_group = anova_table$`Pr(>F)`[1],
    Genotype = anova_table$`Pr(>F)`[2],
    Treatment = anova_table$`Pr(>F)`[3],
    `OD_group*Genotype` = anova_table$`Pr(>F)`[4],
    `OD_group*Treatment` = anova_table$`Pr(>F)`[5],
    `Genotype*Treatment` = anova_table$`Pr(>F)`[6],
    `OD_group*Genotype*Treatment` = anova_table$`Pr(>F)`[7]
  )
  
  # Append the new row to Data_stats_1
  Data_stats_1 <- rbind(Data_stats_1, new_row)
}

Risk_Score_anova <- TableS3[["Risk_Score"]]
Prop_Day4_Alive_anova <- TableS3[["Prop_Day4_Alive"]]
Fs_Day4_anova <- TableS3[["Fs_Day4"]]
PP_Day4_anova <- TableS3[["PP_Day4"]]
Fs_T_anova <- TableS3[["Fs_T"]]
PP_T_anova <- TableS3[["PP_T"]]

Data_stats <- as.data.frame(melt(Data_stats_1, id.vars = c("variable"), variable.name = "vars"))

Data_stats$pval <- ifelse(Data_stats$value>0.05,"Not Significant",
                          ifelse(Data_stats$value<0.05 &
                                   Data_stats$value>0.01, "P<0.05",
                                 ifelse(Data_stats$value<0.01 &
                                          Data_stats$value>0.001, "P<0.01",
                                        ifelse(Data_stats$value<0.001, "P<0.001","NA"))))

custom_labels_vars <- c("OD_group" = "OD", 
                   "Genotype" = "Host genotype", 
                   "Treatment" = "Pathogen species",
                   "OD_group.Genotype" = "OD * Host genotype",
                   "OD_group.Treatment" = "OD * Pathogen species",
                   "Genotype.Treatment" = "Host genotype * Pathogen species",
                   "OD_group.Genotype.Treatment" = "OD * Host genotype * Pathogen species")

custom_labels_variable <- c("Fs_Day4" = "Fs on Day4", 
                   "Fs_T" = "Fs_T", 
                   "PP_Day4" = "PP on Day4",
                   "PP_T" = "PP_T",
                   "Prop_Day4_Alive" = "Proportion alive on Day4",
                   "Risk_Score" = "Risk score")

Fig4A <- ggplot(Data_stats, aes(x=as.factor(vars),
                         y=as.factor(variable),
                         fill=as.factor(pval))) +
  geom_tile() +
  theme_bw() +
  labs(x="Factor",
       y="Measure of pathogenicity", 
       fill="p-value") +
  scale_x_discrete(labels = custom_labels_vars) +
  scale_y_discrete(labels = custom_labels_variable) +
  scale_fill_manual(values = c("white","maroon","red","pink")) +
  theme(axis.text.y = element_text(size=40,
                                   colour = c("darkorchid2","darkorchid2","orange","orange","darkorchid2","darkorchid2")),
        axis.text.x = element_text(size=40, angle = 45, vjust = 1, hjust = 1,
                                   colour = c("black","darkorchid2","orange","black","black","black","black")),
        axis.title = element_text(size = 40),
        legend.position = "top",
        legend.text = element_text(size = 40),
        legend.title = element_text(size = 40))

Figure 4B

# Select the columns for analysis
Corrdata <- Data[,c(27,34,42,46,48,49)]

# Compute the correlation matrix
M <- cor(Corrdata)

TableS4 <- as.data.frame(M)

# Compute the R^2 matrix
R2 <- M^2

# Melt R2 matrix into a long format
df <- melt(R2)

# Filter out the upper triangle of the matrix
df <- df[as.numeric(df$Var1) >= as.numeric(df$Var2), ]

# Correlation plot - Figure 4B
Fig4B <- ggplot(data = df, 
                aes(Var1, Var2, 
                    fill = value)) +
  geom_tile(color = "black") + 
  scale_fill_gradient2(low = "white", 
                       high = "lightblue", 
                       mid = "grey90", 
                       midpoint = 0.5) +
  geom_text(aes(label = round(value, 2)), 
            size = 20) +  
  theme_classic() +  
  theme(axis.text.x = element_text(colour = c("darkorchid2", 
                                              "orange", 
                                              "darkorchid2", 
                                              "darkorchid2", 
                                              "darkorchid2", 
                                              "orange"),
                                   size = 40,
                                   angle = 45, 
                                   vjust = 1, 
                                   hjust = 1,),
        axis.text.y = element_text(colour = c("darkorchid2", 
                                              "orange", 
                                              "darkorchid2", 
                                              "darkorchid2", 
                                              "darkorchid2", 
                                              "orange"),
                                   size = 40),
        legend.position = "none") +
  scale_x_discrete(labels = custom_labels_variable) +  
  scale_y_discrete(labels = custom_labels_variable) +
  labs(x="",
       y="")

Figure4 <- ggarrange(Fig4A, Fig4B, 
                  nrow = 1, ncol = 2, 
                  labels = "AUTO",
                  font.label = list(size = 50))

Figure 4A P-values from ANOVA to indicate the significance of each factor corresponding to each measure of pathogenicity

Figure 4B Heatmap showing R2 (R=correlation coefficient) values to assess the differences between different measures of host susceptibility (survival, risk scores, and fraction symptomatic) and pathogen virulence (pathogen potential)

[Prop_Day4_Alive: Proportion survival on Day 4 post-infection, PP_Day4: Pathogen potential on Day 4, PP_T: Pathogen potential on median survival day, Fs_Day4: Fraction symptomatic on Day 4, Fs_T: Fraction symptomatic on median survival day].

Supplementary Figures

Figure S1

FigureS1 <- ggplot(LoadData, aes(x=OD,
                 y=`Bacterial_Load_CFU/fly`,
                 col=Treatment)) +
  geom_point(alpha=0.5,
             position = position_jitterdodge()) +
  stat_smooth(method = 'lm',
              se=FALSE) +
  theme_bw() +
  theme(axis.text=element_text(size = 24),
        legend.text = element_text(size = 20, face = "italic"),
        legend.position = "top",
        title=element_text(size = 20)) +
  labs(y="Bacterial load (CFU/fly)")

Figure S1- Bacterial load (CFU/fly) vs OD (A600nm) for different ODs and different pathogens. Single infected flies were homogenized in 500microL LB and plated on LB plates immediately after the infection.

Figure S2

#Subsetting Data for PCA
Subset <- Data[,c(1:4,24:37,39:46,48:50)]

row.names(Subset) = paste(Subset$Genotype,Subset$Treatment,Subset$OD_group,Subset$Vial_ID)

PCA <- prcomp(t(Subset[,c(5:29)]), scale. = TRUE)
summary(PCA)

## Importance of components:
##                            PC1     PC2     PC3     PC4       PC5       PC6
## Standard deviation     19.3630 0.24270 0.10800 0.04423 6.701e-07 3.891e-07
## Proportion of Variance  0.9998 0.00016 0.00003 0.00001 0.000e+00 0.000e+00
## Cumulative Proportion   0.9998 0.99996 0.99999 1.00000 1.000e+00 1.000e+00
##                             PC7       PC8       PC9      PC10      PC11
## Standard deviation     3.08e-07 2.709e-07 2.338e-07 1.684e-07 1.573e-07
## Proportion of Variance 0.00e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion  1.00e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
##                             PC12      PC13      PC14      PC15      PC16
## Standard deviation     1.424e-07 8.313e-08 7.632e-08 5.689e-08 4.971e-08
## Proportion of Variance 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
##                             PC17      PC18      PC19      PC20     PC21
## Standard deviation     3.581e-08 2.793e-08 2.373e-08 2.282e-08 1.89e-08
## Proportion of Variance 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.00e+00
## Cumulative Proportion  1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.00e+00
##                             PC22      PC23      PC24      PC25
## Standard deviation     1.308e-08 7.453e-09 5.123e-09 4.838e-16
## Proportion of Variance 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion  1.000e+00 1.000e+00 1.000e+00 1.000e+00

PCA_scores <- PCA$x

PCA_2 <- cbind(Subset,PCA_scores)

FigureS2 <- ggplot(as.data.frame(PCA_scores),
                   aes(x=PC1,
                       y = PC2,
                       col = row.names(PCA_scores))) +
  geom_point(position = position_jitterdodge(),
             size = 3.5, 
             alpha = 0.6) +
  theme_bw() +
  geom_label_repel(aes(label = row.names(PCA_scores)), 
                   nudge_x = -2, 
                   nudge_y = 2,
                   max.overlaps = Inf,
                   box.padding = 0.5,
                   point.padding = 0.3,
                   force = 2,
                   show.legend = FALSE,
                   label.size = 0.5) +
  theme(legend.position = "none",
        axis.text = element_text(size = 30),
        axis.title = element_text(size = 30)) +
  labs(y = "PC2: 23.66%",
       x = "PC1: 62.44%",
       col = "")

Figure S2 Principal component analysis (PCA) over different measures of host susceptibility (survival, risk scores, and fraction symptomatic) and pathogen virulence (pathogen potential). [Prop_DayX_Alive: Proportion survival on Day X post-infection, PP_DayX: Pathogen potential on Day X post-infection, PP_T: Pathogen potential on median survival day, Fs_T/T: Fraction symptomatic on median survival day/ Median survival day, Fs_DayX: Fraction symptomatic on Day X post-infection, Fs_T: Fraction symptomatic on median survival day].

Statistical Analysis

Table S1

Table S1: Statistical analysis for analyzing the effect of OD and treatment (pathogen species) on log-transformed bacterial load. Data was fitted to a linear mixed-effects model.

[Model: lm( log(Bacterial_Load_CFU/fly+1) ~ OD + Treatment)]

Model1 <- lm(`Bacterial_Load_CFU/fly`~OD+Treatment,data = LoadData)

MSum1 <- summary(Model1)
TableS1 <- as.data.frame(MSum1$coefficients)
print(TableS1)

##                      Estimate Std. Error    t value    Pr(>|t|)
## (Intercept)         -173.1791  1319.6811 -0.1312280 0.896057039
## OD                   976.3464   281.7327  3.4655067 0.001012516
## TreatmentP. sneebia  191.8333  1632.8846  0.1174813 0.906891692

Table S2

Table S2: Statistical analysis for analyzing the effect of OD, treatment (pathogen species), and host genotype on survival. Data was fitted to a mixed-effects cox-proportional hazards model.

[Model: coxph(Survival ~ OD * Treatment * Genotype]

cox <- coxph(Surv(Day_Dead, censor) ~ OD_group * Treatment * Genotype +
               frailty(Date),
             data = Data_1)

cox_summary <- summary(cox)
TableS2 <- as.data.frame(cox_summary$coefficients)


print(TableS2)

##                                      coef     se(coef)          se2
## OD_group0.01                  2.465579129 4.924019e-01 4.882623e-01
## OD_group0.1                   2.543120127 4.850304e-01 4.808286e-01
## OD_group0.5                   2.715388454 4.744765e-01 4.744764e-01
## OD_group1                    -0.141432999 4.765642e-01 4.765642e-01
## OD_group4                     2.853860034 4.698408e-01 4.698407e-01
## OD_group8                     4.149345144 4.736175e-01 4.693276e-01
## TreatmentP..alcalifaciens    -1.671417849 6.982090e+05 6.982090e+05
## TreatmentP..burhodogranar    -1.914697543 2.427268e-01 2.328724e-01
## TreatmentP..rettgeri         -1.246073645 3.722759e+05 3.722759e+05
## TreatmentP..sneebia                    NA 0.000000e+00 0.000000e+00
## GenotypeDptR                 -0.382912516 6.732125e-01 6.712090e-01
## GenotypeDptS                 -0.309707667 6.717803e-01 6.709762e-01
## frailty.Date.                          NA           NA           NA
## OD_group0.01.TreatmentP..     2.399820178 6.982090e+05 6.982090e+05
## OD_group0.1.TreatmentP..a     2.156073003 6.982090e+05 6.982090e+05
## OD_group0.5.TreatmentP..a     2.165868733 6.982090e+05 6.982090e+05
## OD_group1.TreatmentP..alc     5.073849259 6.982090e+05 6.982090e+05
## OD_group4.TreatmentP..alc     2.598092750 6.982090e+05 6.982090e+05
## OD_group8.TreatmentP..alc     1.350307927 6.982090e+05 6.982090e+05
## OD_group0.01.TreatmentP...1   1.371817724 3.648377e-01 3.648376e-01
## OD_group0.1.TreatmentP..b     0.752297228 3.760251e-01 3.760251e-01
## OD_group0.5.TreatmentP..b     0.541166546 3.772288e-01 3.718072e-01
## OD_group1.TreatmentP..bur     4.213991301 3.482430e-01 3.423709e-01
## OD_group4.TreatmentP..bur     0.996788604 3.425627e-01 3.365727e-01
## OD_group8.TreatmentP..bur              NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...2   0.562711101 3.722759e+05 3.722759e+05
## OD_group0.1.TreatmentP..r     1.452804529 3.722759e+05 3.722759e+05
## OD_group0.5.TreatmentP..r     1.502594030 3.722759e+05 3.722759e+05
## OD_group1.TreatmentP..ret     4.785627270 3.722759e+05 3.722759e+05
## OD_group4.TreatmentP..ret     1.264499761 3.722759e+05 3.722759e+05
## OD_group8.TreatmentP..ret     0.171801966 3.722759e+05 3.722759e+05
## OD_group0.01.TreatmentP...3            NA 0.000000e+00 0.000000e+00
## OD_group0.1.TreatmentP..s              NA 0.000000e+00 0.000000e+00
## OD_group0.5.TreatmentP..s              NA 0.000000e+00 0.000000e+00
## OD_group1.TreatmentP..sne              NA 0.000000e+00 0.000000e+00
## OD_group4.TreatmentP..sne              NA 0.000000e+00 0.000000e+00
## OD_group8.TreatmentP..sne              NA 0.000000e+00 0.000000e+00
## OD_group0.01.GenotypeDptR     0.470713354 7.281076e-01 7.262556e-01
## OD_group0.1.GenotypeDptR      0.883859006 7.142435e-01 7.123480e-01
## OD_group0.5.GenotypeDptR      1.055857403 7.045312e-01 7.045311e-01
## OD_group1.GenotypeDptR        0.937833652 7.075987e-01 7.070637e-01
## OD_group4.GenotypeDptR        0.452329473 7.009321e-01 7.003985e-01
## OD_group8.GenotypeDptR        0.420449012 7.025855e-01 7.006654e-01
## OD_group0.01.GenotypeDptS     0.449644062 7.240409e-01 7.232943e-01
## OD_group0.1.GenotypeDptS      0.266941237 7.143786e-01 7.136220e-01
## OD_group0.5.GenotypeDptS      0.882105746 7.048471e-01 7.048330e-01
## OD_group1.GenotypeDptS        1.035919516 7.059839e-01 7.059634e-01
## OD_group4.GenotypeDptS        0.121836449 7.015740e-01 7.015551e-01
## OD_group8.GenotypeDptS       -0.060228758 7.018380e-01 7.010733e-01
## TreatmentP..alcalifaciens.1  -0.085236782 2.835734e-01 2.835734e-01
## TreatmentP..burhodogranar.1  -0.042165244 3.289776e-01 3.289776e-01
## TreatmentP..rettgeri.Geno     0.004937081 2.843097e-01 2.843097e-01
## TreatmentP..sneebia.Genot              NA 0.000000e+00 0.000000e+00
## TreatmentP..alcalifaciens.2   0.322236139 2.851122e-01 2.851122e-01
## TreatmentP..burhodogranar.2   0.276863162 3.338959e-01 3.338959e-01
## TreatmentP..rettgeri.Geno.1  -0.187161768 2.937825e-01 2.937825e-01
## TreatmentP..sneebia.Genot.1            NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...4  -0.314937199 4.504374e-01 4.504374e-01
## OD_group0.1.TreatmentP..a.1  -0.330612917 4.257572e-01 4.257572e-01
## OD_group0.5.TreatmentP..a.1  -0.106091002 4.127077e-01 4.094277e-01
## OD_group1.TreatmentP..alc.1   0.097547468 4.145128e-01 4.137762e-01
## OD_group4.TreatmentP..alc.1  -0.029708361 4.019184e-01 4.011670e-01
## OD_group8.TreatmentP..alc.1            NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...5  -0.710053077 5.363291e-01 5.363291e-01
## OD_group0.1.TreatmentP..b.1   0.158952671 5.076546e-01 5.076546e-01
## OD_group0.5.TreatmentP..b.1  -0.847586508 5.366671e-01 5.341718e-01
## OD_group1.TreatmentP..bur.1  -2.066817617 5.528053e-01 5.522705e-01
## OD_group4.TreatmentP..bur.1   0.259844109 4.714907e-01 4.708480e-01
## OD_group8.TreatmentP..bur.1            NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...6   0.224937309 4.867779e-01 4.867779e-01
## OD_group0.1.TreatmentP..r.1  -0.304239563 4.286863e-01 4.286863e-01
## OD_group0.5.TreatmentP..r.1  -0.565978533 4.137873e-01 4.105542e-01
## OD_group1.TreatmentP..ret.1  -0.677703633 4.150179e-01 4.143087e-01
## OD_group4.TreatmentP..ret.1  -0.272944117 4.111758e-01 4.104495e-01
## OD_group8.TreatmentP..ret.1            NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...7            NA 0.000000e+00 0.000000e+00
## OD_group0.1.TreatmentP..s.1            NA 0.000000e+00 0.000000e+00
## OD_group0.5.TreatmentP..s.1            NA 0.000000e+00 0.000000e+00
## OD_group1.TreatmentP..sne.1            NA 0.000000e+00 0.000000e+00
## OD_group4.TreatmentP..sne.1            NA 0.000000e+00 0.000000e+00
## OD_group8.TreatmentP..sne.1            NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...8  -1.431941830 4.582501e-01 4.525892e-01
## OD_group0.1.TreatmentP..a.2  -0.537402197 4.368553e-01 4.309135e-01
## OD_group0.5.TreatmentP..a.2  -0.826002197 4.170152e-01 4.098035e-01
## OD_group1.TreatmentP..alc.2  -1.058159287 4.204203e-01 4.133327e-01
## OD_group4.TreatmentP..alc.2  -0.265166412 4.054168e-01 4.044784e-01
## OD_group8.TreatmentP..alc.2            NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...9  -1.533479745 5.548007e-01 5.548007e-01
## OD_group0.1.TreatmentP..b.2  -0.542555529 5.454616e-01 5.454616e-01
## OD_group0.5.TreatmentP..b.2  -1.013040316 5.336373e-01 5.328856e-01
## OD_group1.TreatmentP..bur.2  -2.096057726 5.177448e-01 5.170239e-01
## OD_group4.TreatmentP..bur.2   0.043710955 4.796295e-01 4.788321e-01
## OD_group8.TreatmentP..bur.2            NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...10  0.030270044 4.943688e-01 4.943688e-01
## OD_group0.1.TreatmentP..r.2  -0.111723849 4.428810e-01 4.428810e-01
## OD_group0.5.TreatmentP..r.2  -0.646098913 4.222030e-01 4.212477e-01
## OD_group1.TreatmentP..ret.2  -0.975086453 4.206219e-01 4.197294e-01
## OD_group4.TreatmentP..ret.2  -0.366886886 4.288274e-01 4.279366e-01
## OD_group8.TreatmentP..ret.2            NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...11           NA 0.000000e+00 0.000000e+00
## OD_group0.1.TreatmentP..s.2            NA 0.000000e+00 0.000000e+00
## OD_group0.5.TreatmentP..s.2            NA 0.000000e+00 0.000000e+00
## OD_group1.TreatmentP..sne.2            NA 0.000000e+00 0.000000e+00
## OD_group4.TreatmentP..sne.2            NA 0.000000e+00 0.000000e+00
## OD_group8.TreatmentP..sne.2            NA 0.000000e+00 0.000000e+00
##                                     Chisq          DF            p
## OD_group0.01                 2.507255e+01  1.00000000 5.521322e-07
## OD_group0.1                  2.749134e+01  1.00000000 1.577995e-07
## OD_group0.5                  3.275175e+01  1.00000000 1.047119e-08
## OD_group1                    8.807622e-02  1.00000000 7.666372e-01
## OD_group4                    3.689469e+01  1.00000000 1.246852e-09
## OD_group8                    7.675448e+01  1.00000000 1.935887e-18
## TreatmentP..alcalifaciens    5.730588e-12  1.00000000 9.999981e-01
## TreatmentP..burhodogranar    6.222500e+01  1.00000000 3.063741e-15
## TreatmentP..rettgeri         1.120360e-11  1.00000000 9.999973e-01
## TreatmentP..sneebia                    NA  1.00000000           NA
## GenotypeDptR                 3.235153e-01  1.00000000 5.695027e-01
## GenotypeDptS                 2.125442e-01  1.00000000 6.447805e-01
## frailty.Date.                3.525123e+00 -0.08584406 2.353071e-02
## OD_group0.01.TreatmentP..    1.181372e-11  1.00000000 9.999973e-01
## OD_group0.1.TreatmentP..a    9.535776e-12  1.00000000 9.999975e-01
## OD_group0.5.TreatmentP..a    9.622621e-12  1.00000000 9.999975e-01
## OD_group1.TreatmentP..alc    5.280855e-11  1.00000000 9.999942e-01
## OD_group4.TreatmentP..alc    1.384645e-11  1.00000000 9.999970e-01
## OD_group8.TreatmentP..alc    3.740199e-12  1.00000000 9.999985e-01
## OD_group0.01.TreatmentP...1  1.413818e+01  1.00000000 1.698613e-04
## OD_group0.1.TreatmentP..b    4.002628e+00  1.00000000 4.542937e-02
## OD_group0.5.TreatmentP..b    2.058032e+00  1.00000000 1.514056e-01
## OD_group1.TreatmentP..bur    1.464275e+02  1.00000000 1.046858e-33
## OD_group4.TreatmentP..bur    8.466930e+00  1.00000000 3.616612e-03
## OD_group8.TreatmentP..bur              NA  1.00000000           NA
## OD_group0.01.TreatmentP...2  2.284763e-12  1.00000000 9.999988e-01
## OD_group0.1.TreatmentP..r    1.522946e-11  1.00000000 9.999969e-01
## OD_group0.5.TreatmentP..r    1.629121e-11  1.00000000 9.999968e-01
## OD_group1.TreatmentP..ret    1.652524e-10  1.00000000 9.999897e-01
## OD_group4.TreatmentP..ret    1.153739e-11  1.00000000 9.999973e-01
## OD_group8.TreatmentP..ret    2.129739e-13  1.00000000 9.999996e-01
## OD_group0.01.TreatmentP...3            NA  1.00000000           NA
## OD_group0.1.TreatmentP..s              NA  1.00000000           NA
## OD_group0.5.TreatmentP..s              NA  1.00000000           NA
## OD_group1.TreatmentP..sne              NA  1.00000000           NA
## OD_group4.TreatmentP..sne              NA  1.00000000           NA
## OD_group8.TreatmentP..sne              NA  1.00000000           NA
## OD_group0.01.GenotypeDptR    4.179477e-01  1.00000000 5.179629e-01
## OD_group0.1.GenotypeDptR     1.531346e+00  1.00000000 2.159105e-01
## OD_group0.5.GenotypeDptR     2.246002e+00  1.00000000 1.339601e-01
## OD_group1.GenotypeDptR       1.756619e+00  1.00000000 1.850468e-01
## OD_group4.GenotypeDptR       4.164453e-01  1.00000000 5.187161e-01
## OD_group8.GenotypeDptR       3.581198e-01  1.00000000 5.495523e-01
## OD_group0.01.GenotypeDptS    3.856661e-01  1.00000000 5.345863e-01
## OD_group0.1.GenotypeDptS     1.396286e-01  1.00000000 7.086505e-01
## OD_group0.5.GenotypeDptS     1.566215e+00  1.00000000 2.107575e-01
## OD_group1.GenotypeDptS       2.153091e+00  1.00000000 1.422831e-01
## OD_group4.GenotypeDptS       3.015835e-02  1.00000000 8.621314e-01
## OD_group8.GenotypeDptS       7.364344e-03  1.00000000 9.316129e-01
## TreatmentP..alcalifaciens.1  9.034895e-02  1.00000000 7.637340e-01
## TreatmentP..burhodogranar.1  1.642770e-02  1.00000000 8.980140e-01
## TreatmentP..rettgeri.Geno    3.015486e-04  1.00000000 9.861453e-01
## TreatmentP..sneebia.Genot              NA  1.00000000           NA
## TreatmentP..alcalifaciens.2  1.277370e+00  1.00000000 2.583886e-01
## TreatmentP..burhodogranar.2  6.875560e-01  1.00000000 4.069969e-01
## TreatmentP..rettgeri.Geno.1  4.058658e-01  1.00000000 5.240754e-01
## TreatmentP..sneebia.Genot.1            NA  1.00000000           NA
## OD_group0.01.TreatmentP...4  4.888538e-01  1.00000000 4.844391e-01
## OD_group0.1.TreatmentP..a.1  6.029975e-01  1.00000000 4.374366e-01
## OD_group0.5.TreatmentP..a.1  6.608028e-02  1.00000000 7.971318e-01
## OD_group1.TreatmentP..alc.1  5.538041e-02  1.00000000 8.139523e-01
## OD_group4.TreatmentP..alc.1  5.463633e-03  1.00000000 9.410769e-01
## OD_group8.TreatmentP..alc.1            NA  1.00000000           NA
## OD_group0.01.TreatmentP...5  1.752746e+00  1.00000000 1.855319e-01
## OD_group0.1.TreatmentP..b.1  9.803904e-02  1.00000000 7.541957e-01
## OD_group0.5.TreatmentP..b.1  2.494354e+00  1.00000000 1.142552e-01
## OD_group1.TreatmentP..bur.1  1.397848e+01  1.00000000 1.849152e-04
## OD_group4.TreatmentP..bur.1  3.037243e-01  1.00000000 5.815570e-01
## OD_group8.TreatmentP..bur.1            NA  1.00000000           NA
## OD_group0.01.TreatmentP...6  2.135312e-01  1.00000000 6.440137e-01
## OD_group0.1.TreatmentP..r.1  5.036770e-01  1.00000000 4.778889e-01
## OD_group0.5.TreatmentP..r.1  1.870878e+00  1.00000000 1.713745e-01
## OD_group1.TreatmentP..ret.1  2.666527e+00  1.00000000 1.024794e-01
## OD_group4.TreatmentP..ret.1  4.406487e-01  1.00000000 5.068095e-01
## OD_group8.TreatmentP..ret.1            NA  1.00000000           NA
## OD_group0.01.TreatmentP...7            NA  1.00000000           NA
## OD_group0.1.TreatmentP..s.1            NA  1.00000000           NA
## OD_group0.5.TreatmentP..s.1            NA  1.00000000           NA
## OD_group1.TreatmentP..sne.1            NA  1.00000000           NA
## OD_group4.TreatmentP..sne.1            NA  1.00000000           NA
## OD_group8.TreatmentP..sne.1            NA  1.00000000           NA
## OD_group0.01.TreatmentP...8  9.764399e+00  1.00000000 1.779236e-03
## OD_group0.1.TreatmentP..a.2  1.513295e+00  1.00000000 2.186370e-01
## OD_group0.5.TreatmentP..a.2  3.923364e+00  1.00000000 4.761949e-02
## OD_group1.TreatmentP..alc.2  6.334828e+00  1.00000000 1.183896e-02
## OD_group4.TreatmentP..alc.2  4.277930e-01  1.00000000 5.130739e-01
## OD_group8.TreatmentP..alc.2            NA  1.00000000           NA
## OD_group0.01.TreatmentP...9  7.639801e+00  1.00000000 5.709422e-03
## OD_group0.1.TreatmentP..b.2  9.893728e-01  1.00000000 3.198957e-01
## OD_group0.5.TreatmentP..b.2  3.603803e+00  1.00000000 5.764755e-02
## OD_group1.TreatmentP..bur.2  1.638985e+01  1.00000000 5.156051e-05
## OD_group4.TreatmentP..bur.2  8.305557e-03  1.00000000 9.273855e-01
## OD_group8.TreatmentP..bur.2            NA  1.00000000           NA
## OD_group0.01.TreatmentP...10 3.749074e-03  1.00000000 9.511763e-01
## OD_group0.1.TreatmentP..r.2  6.363818e-02  1.00000000 8.008354e-01
## OD_group0.5.TreatmentP..r.2  2.341831e+00  1.00000000 1.259415e-01
## OD_group1.TreatmentP..ret.2  5.374061e+00  1.00000000 2.043834e-02
## OD_group4.TreatmentP..ret.2  7.319801e-01  1.00000000 3.922421e-01
## OD_group8.TreatmentP..ret.2            NA  1.00000000           NA
## OD_group0.01.TreatmentP...11           NA  1.00000000           NA
## OD_group0.1.TreatmentP..s.2            NA  1.00000000           NA
## OD_group0.5.TreatmentP..s.2            NA  1.00000000           NA
## OD_group1.TreatmentP..sne.2            NA  1.00000000           NA
## OD_group4.TreatmentP..sne.2            NA  1.00000000           NA
## OD_group8.TreatmentP..sne.2            NA  1.00000000           NA

Table S3

Table S3: ANOVA for analyzing the effect of OD, treatment (pathogen species), and host genotype on different measures of pathogenicity, corresponding to Figure 4A.

[Model: aov(“Measure of pathogenicity” ~ OD * Treatment * Genotype]

# Variables of interest
variables_of_interest <- c("Risk_Score", 
                           "Prop_Day4_Alive", 
                           "Fs_Day4", 
                           "PP_Day4", 
                           "Fs_T", 
                           "PP_T")

# Create an empty dataframe to store p-values matrices
Data_stats_1 <- data.frame(
  variable = character(),
  OD_group = character(),
  Genotype = character(),
  Treatment = character(),
  `OD_group*Genotype` = character(),
  `OD_group*Treatment` = character(),
  `Genotype*Treatment` = character(),
  `OD_group*Genotype*Treatment` = character(),
  stringsAsFactors = FALSE
)
TableS3 <- list()

# Loop through each variable of interest
for (variable in variables_of_interest) {
  # Fit an ANOVA model
  model <- aov(as.formula(paste(variable, "~ OD_group * Genotype * Treatment + Error(Vial_ID+Date)")), data = Data)
  
  # Extract ANOVA table
  anova_table <- as.data.frame(summary(model$Within)[[1]])
  
  TableS3[[variable]] <- anova_table
  
  # Create a new row to append
  new_row <- data.frame(
    variable = variable,
    OD_group = anova_table$`Pr(>F)`[1],
    Genotype = anova_table$`Pr(>F)`[2],
    Treatment = anova_table$`Pr(>F)`[3],
    `OD_group*Genotype` = anova_table$`Pr(>F)`[4],
    `OD_group*Treatment` = anova_table$`Pr(>F)`[5],
    `Genotype*Treatment` = anova_table$`Pr(>F)`[6],
    `OD_group*Genotype*Treatment` = anova_table$`Pr(>F)`[7]
  )
  
  # Append the new row to Data_stats_1
  Data_stats_1 <- rbind(Data_stats_1, new_row)
}


Risk_Score_anova <- TableS3[["Risk_Score"]]
Prop_Day4_Alive_anova <- TableS3[["Prop_Day4_Alive"]]
Fs_Day4_anova <- TableS3[["Fs_Day4"]]
PP_Day4_anova <- TableS3[["PP_Day4"]]
Fs_T_anova <- TableS3[["Fs_T"]]
PP_T_anova <- TableS3[["PP_T"]]

A. Risk score

##                              Df       Sum Sq      Mean Sq     F value
## OD_group                      1 1.721013e+17 1.721013e+17 188.7850978
## Genotype                      2 2.380668e+16 1.190334e+16  13.0572697
## Treatment                     1 8.631467e+16 8.631467e+16  94.6821660
## OD_group:Genotype             2 8.151767e+14 4.075883e+14   0.4471007
## OD_group:Treatment            3 3.397697e+16 1.132566e+16  12.4235869
## Genotype:Treatment            8 2.442590e+16 3.053237e+15   3.3492234
## OD_group:Genotype:Treatment   6 2.156488e+16 3.594146e+15   3.9425694
## Residuals                   343 3.126875e+17 9.116254e+14          NA
##                                   Pr(>F)
## OD_group                    1.569612e-34
## Genotype                    3.426131e-06
## Treatment                   6.399903e-20
## OD_group:Genotype           6.398517e-01
## OD_group:Treatment          9.870310e-08
## Genotype:Treatment          1.041046e-03
## OD_group:Genotype:Treatment 7.938854e-04
## Residuals                             NA

B. Proportion alive on day 4

##                              Df     Sum Sq    Mean Sq     F value       Pr(>F)
## OD_group                      1 2.84511946 2.84511946 146.4453151 2.580869e-28
## Genotype                      2 0.22638178 0.11319089   5.8262143 3.248910e-03
## Treatment                     1 6.24734970 6.24734970 321.5664960 3.371913e-51
## OD_group:Genotype             2 0.00645126 0.00322563   0.1660311 8.470879e-01
## OD_group:Treatment            3 0.36332099 0.12110700   6.2336758 3.933111e-04
## Genotype:Treatment            8 0.27880303 0.03485038   1.7938350 7.717203e-02
## OD_group:Genotype:Treatment   6 0.17242848 0.02873808   1.4792198 1.843522e-01
## Residuals                   343 6.66375687 0.01942786          NA           NA

C. Pathogen potential on median survival day (PP_T)

##                              Df      Sum Sq      Mean Sq      F value
## OD_group                      1 192.2315861 192.23158606 119.93771288
## Genotype                      2   0.7081393   0.35406963   0.22091219
## Treatment                     1  27.7883581  27.78835812  17.33779649
## OD_group:Genotype             2   0.1805177   0.09025885   0.05631457
## OD_group:Treatment            3  34.4844024  11.49480081   7.17187091
## Genotype:Treatment            8  15.3005495   1.91256869   1.19329565
## OD_group:Genotype:Treatment   6  12.7716303   2.12860506   1.32808571
## Residuals                   343 549.7473016   1.60276181           NA
##                                   Pr(>F)
## OD_group                    3.886453e-24
## Genotype                    8.019011e-01
## Treatment                   3.962438e-05
## OD_group:Genotype           9.452505e-01
## OD_group:Treatment          1.105386e-04
## Genotype:Treatment          3.020478e-01
## OD_group:Genotype:Treatment 2.437902e-01
## Residuals                             NA

D. Pathogen potential on day 4 (PP_Day4)

##                              Df      Sum Sq      Mean Sq      F value
## OD_group                      1 419.1892066 419.18920656 236.74644014
## Genotype                      2   1.7895488   0.89477438   0.50534376
## Treatment                     1  72.7881222  72.78812220  41.10871307
## OD_group:Genotype             2   0.1078957   0.05394785   0.03046825
## OD_group:Treatment            3  68.1946760  22.73155868  12.83815402
## Genotype:Treatment            8  27.7610991   3.47013739   1.95983737
## OD_group:Genotype:Treatment   6  28.9251722   4.82086203   2.72268920
## Residuals                   343 607.3244344   1.77062517           NA
##                                   Pr(>F)
## OD_group                    5.423226e-41
## Genotype                    6.037466e-01
## Treatment                   4.768189e-10
## OD_group:Genotype           9.699939e-01
## OD_group:Treatment          5.717453e-08
## Genotype:Treatment          5.072378e-02
## OD_group:Genotype:Treatment 1.347991e-02
## Residuals                             NA

E. Fraction symptomatic on median survival day (Fs_T)

##                              Df      Sum Sq     Mean Sq    F value       Pr(>F)
## OD_group                      1 0.208976905 0.208976905 11.0224812 9.968983e-04
## Genotype                      2 0.117802824 0.058901412  3.1067534 4.600759e-02
## Treatment                     1 0.327587015 0.327587015 17.2785683 4.081741e-05
## OD_group:Genotype             2 0.000835166 0.000417583  0.0220254 9.782168e-01
## OD_group:Treatment            3 0.163925138 0.054641713  2.8820757 3.588904e-02
## Genotype:Treatment            8 0.354809753 0.044351219  2.3393039 1.855107e-02
## OD_group:Genotype:Treatment   6 0.184355977 0.030725996  1.6206418 1.404555e-01
## Residuals                   343 6.502989371 0.018959153         NA           NA

F. Fraction symptomatic on day 4 (Fs_Day4)

##                              Df     Sum Sq     Mean Sq    F value       Pr(>F)
## OD_group                      1 0.39252080 0.392520798 29.3149653 1.158839e-07
## Genotype                      2 0.05260709 0.026303544  1.9644500 1.418076e-01
## Treatment                     1 0.98637215 0.986372151 73.6660720 3.264645e-16
## OD_group:Genotype             2 0.01769875 0.008849373  0.6609053 5.170400e-01
## OD_group:Treatment            3 0.05409545 0.018031818  1.3466856 2.590671e-01
## Genotype:Treatment            8 0.24379423 0.030474279  2.2759366 2.202909e-02
## OD_group:Genotype:Treatment   6 0.16760849 0.027934748  2.0862746 5.428488e-02
## Residuals                   343 4.59269293 0.013389775         NA           NA

Table S4

Table S4: Correlation matrix for different measures of pathogenicity containing R values corresponding to Figure 4B.

# Select the columns for analysis
Corrdata <- Data[,c(27,34,42,46,48:50)]

# Compute the correlation matrix
M <- cor(Corrdata)

TableS4 <- as.data.frame(M)

print(TableS4)

##                    Fs_Day4    PP_Day4 Prop_Day4_Alive Risk_Score       Fs_T
## Fs_Day4          1.0000000  0.3035784      -0.7881867  0.6467606  0.6756411
## PP_Day4          0.3035784  1.0000000      -0.3719103  0.2167741  0.2059444
## Prop_Day4_Alive -0.7881867 -0.3719103       1.0000000 -0.8139554 -0.5490079
## Risk_Score       0.6467606  0.2167741      -0.8139554  1.0000000  0.6024717
## Fs_T             0.6756411  0.2059444      -0.5490079  0.6024717  1.0000000
## PP_T             0.1941431  0.8643599      -0.2450816  0.1913188  0.3243466
## Fs_T/T           0.6789812  0.3029858      -0.8681956  0.8430373  0.5795366
##                       PP_T     Fs_T/T
## Fs_Day4          0.1941431  0.6789812
## PP_Day4          0.8643599  0.3029858
## Prop_Day4_Alive -0.2450816 -0.8681956
## Risk_Score       0.1913188  0.8430373
## Fs_T             0.3243466  0.5795366
## PP_T             1.0000000  0.2066825
## Fs_T/T           0.2066825  1.0000000